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In this paper we present a perturbative procedure that allows one to numerically solve diffusive 
non-Markovian Stochastic Schrodinger equations, for a wide range of memory functions. To illustrate 
this procedure numerical results are presented for a classically driven two level atom immersed in a 
environment with a simple memory function. It is observed that as the order of the perturbation is 
increased the numerical results for the ensembled average state p re d(t) approach the exact reduced 
state found via Imamoglu's enlarged system method [Phys. Rev. A. 50, 3650 (1994)]. 
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I. INTRODUCTION 

A common problem in physics is to model open quan- 
tum systems. They consists of a small system immersed 
in a bath (environment). Due to the large Hilbert space 
of the bath it is convenient to describe the system by its 
reduced state. The reduced state is defined as 



p led (t) =Ti- bath [|#(*))(tf(t)|]. 



(1.1) 



where \^(t)) is the combined system state, found from 
the Schrodinger equation for the open quantum system. 

It has been shown [1, 2] by a projection-operator 
method that we can write a general master equation for 
the reduced state as 



-~[H(t),p Ied (t)\ 



/C(i, s)[L]p Ted (s)ds 



(1.2) 

where IC(t, s)[L] is the 'memory time' superoperator. It 
(operators on) the system operator L and represents how 
the bath affects the system. The problem with this equa- 
tion is that in general IC(t, s)[L] can not be explicitly 
evaluated. 

The most notable approximation used is the Born- 
Markov one. This arises when the environmental influ- 
ences on the system are instantaneous. Mathematical 
consistency requires that this results in a Lindblad mas- 
ter equation, of the form [3] 



For the non-Markovian situation there have been many 
attempts at finding solutions to Eq. (1.2). However, some 
have the problem that it is hard to ensure the positivity 
requirement on p Tcd (t) [8]. A method that does ensure 
the positivity requirement on the reduced state is the 
non-Markovian stochastic Schrodinger equation (SSE) 
approach [9, 10, 11, 12, 13, 14, 15, 16]. A non-Markovian 
SSE generates stochastic pure states \ip z (t)} that should 
satisfy 



Pred (t) = E[\Mt))(Mt)\h 



(1.5) 



Pred(i) 



-~[H(t),p Ied (t)]+jV[L]p Ted (t), 



(1.3) 



where z(t) is some noise function which is non- white and 
E denotes the ensemble average over z(t). To solve these 
non-Markovian SSE one has to take into account the past 
behavior of the system and bath, giving rise to a func- 
tional derivative in the attempt to derive a SSE. This 
presents a problem as for most systems an exact solu- 
tion to the functional derivative does not exists. Thus at 
present an exact non-Markovian SSE only exists for sim- 
ple systems, which can be solve exactly via other meth- 
ods, like the undriven two level atom (TLA). For this and 
more examples see Rcf. [11, 16]. 

Recently Yu, Diosi, Gisin and Strunz (YDGS) have 
developed explicitly a 'post-Markovian' perturbation 
method to first order that allows solutions for systems 
that are close to the Markovian limit [17, 18]. In this pa- 
per we present a perturbation method that can be carried 
to arbitrary order and so is not limited to the post Marko- 
vian regime. However we must place a requirement on 
the form of the memory functions. This requirement is 
that the memory function must take the form 



where T>[L] is the superoperator that represent the damp- 
ing of the system into the bath. It has the form 

T>[L]p Icd = Lp Icd L f - \DLp rcd - \p rcd DL. (1.4) 

This equation can be solved deterministically [4] or by 
the stochastic Schrodinger approach [4, 5, 6, 7]. 



a(t - s) = V |G 



(1.6) 
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for some finite (and, in practice, relatively small) J. It 
should be noted also that we have not proven convergence 
of our perturbation theory and this theory is only valid 
for a zero-temperature bath. 

The format of this paper is as follows. In Sec. II we 
present a general outline of the theory of non-Markovian 
SSE. This is basically a summary of the results of Rcfs. 
[9, 10, 11, 12, 16]. In Sec. Ill our perturbation method 
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is presented. In Sees. IV we outline Imamoglu enlarged 
system method [19, 20]. In Sec. V we apply our per- 
turbation method to a simple system, a driven TLA and 
compare our results for p TC d(t) with the enlarged sys- 
tem methods. In Sec. VI we investigate YDGS post- 
Markovian perturbation method [17, 18]. Finally we con- 
clude with a discussion of the potential applications of 
our results in Sec. VII. 



II. NON-MARKOVIAN STOCHASTIC 
SCHRODINGER EQUATIONS 

In this section we will present an outline of the theory 
we presented in [16], which is an extension of Diosi, Gisin 
and Strunz (DGS) diffusive Non-Markovian SSEs [9, 10, 
11, 12] which allows for real- valued noise z(t). 

A. Underlying Dynamics 

The non-Markovian SSEs developed in references [9, 
10, 11, 12, 16] are valid when the dynamics of the open 
quantum system can be described by the total Hamilto- 
nian 

fftot = H sys ® i + i ® H hat h + V. (2.1) 

The system Hamiltonian is H sys = Hq + H. The bath is 
modelled by a collection of harmonic oscillators, so the 
Hamiltonian for the bath is 

H ha ,th = ^y^^fcajafc, (2.2) 

k 

where and uok are the lowering operator and angular 
frequency of the k th mode respectively. This is the stan- 
dard model for the electromagnetic field. The interaction 
Hamiltonian has the form 

V = ih(W - Db), (2.3) 

where we have defined the bath lowering operators b as 
b = gk&k- That is, the coupling amplitude of the fc th 
mode to the system is g^. 

For calculation purposes we define the non-Markovian 
SSE in an interaction picture. This allows us to move the 
fast dynamics placed on the state by the Hamiltonians 
Hq and -ffbath to the operators. The unitary evolution 
operator for this transformations is 

U(t,0) — e ~K-(^®i+i® i ^b'«th)(t-o)^ ^.4) 

Thus the combined state in the interaction picture is de- 
fine as 

|*(t)) = Crt(t,0)|*(t)sch), (2.5) 
and an arbitrary operator A becomes 

A int = U^(t,0)lU(t,0). (2.6) 



This allows us to write the Schrodingcr equation as 

*!*(*)> = -£[#tat(t) + VU(t)]\*(t)), (2.7) 
where the Hamiltonians are 

H int (t) = tf(t,0)HU(t,0), (2.8) 

and 

V int (t) = ih[Le- iQ % t (t) - Lte^SfctC*)], (2.9) 

where 

knt(t) = J29ka k e- iUkt . (2.10) 

k 

Here we have finally restricted the form of Hq to be such 
that L in the interaction picture simply rotates in the 
complex plane at frequency f2. That is L- mi {t) = Le~ int . 

B. Non-Markovian SSE-Defined 

A non-Markovian SSE is a stochastic differential equa- 
tion for the system state vector \ip z (t)) containing some 
non- white noise z(t). It has the property that when 
\4>z(t))(^Jz(t)\ is averaged over all possible z(t) one ob- 
tains p r ed(t). It should be noted that for a single p rc d(t), 
z(t) can take many different functional forms, and we 
label these different forms as stochastic unravelings [16]. 

In Ref. [16] we showed that non-Markovian SSEs can 
be derived from quantum measurement theory (QMT), 
where the different unravelings correspond to different 
measurements on the bath. The two unravelings we con- 
sidered were the 'coherent' or 'DGS' [9, 10, 11, 12] un- 
raveling and the 'quadrature' unraveling. A special case 
of our quadrature unraveling was published in Ref. [21]. 

As in the Markov limit we can define (at least) two non- 
Markovian SSEs, for each unraveling: one for z(t) chosen 
from an ostensible distribution (a guessed distribution) 
and the other for its actual distribution. The former 
gives a non-Markovian SSE linear in the unnormalised 
state \ip z (t)), while the latter gives a non-Markovian SSE 
non-linear in the normalized state \ip z (t)). In Ref. [16] 
we came to the conclusion that the solution of the actual 
non-Markovian SSE at time t gives the state the system 
will be in if a measurement of the bath is performed at 
that time. Unlike in the Markov case, linking of the 
states through time to make a trajectory turns out to be 
a convenient fiction. However, it has been suggested that 
such trajectories can be given an interpretation within a 
non-standard QMT [22, 23]. 

1. Coherent Unravelling- Outlined 

The first unravelling we consider is the 'coherent' un- 
ravelling. This unravelling arises when the bath is pro- 
jected into a coherent state. We define the coherent state 
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as 



iK})=n^" K|2/2 E^=i 



\n k ), (2.11) 



so that i = l\k I \a k )(a k \d 2 a k . 

In a measurement we can define an operator for the 
measurement process, the noise operator. For this mea- 
surement it must have the coherent basis as its eigenstate, 
so the noise operator is 



z(t) = b iat (t)e 



tin 



(2.12) 



where Q k — — ^- The noise function (eigenvalue of 
the noise operator) is 



(2.13) 



where a k are the results of the projection in the coherent 
basis. 

If we assume an ostensible distribution for a k as being 
the overlap of the coherent state with vacuum state, that 
is, it has the form 

A({a k }) = ({0 k }\{a k })({a k }\{0 k }) = n- K e-^^ 2 , 

(2.14) 

where K = J2k- With this ostensible distribution the 
noise function has the following correlations 

E[z{t)z*(s)] = a(t-s), (2.15a) 
E[z{t)z(s)} = 0. (2.15b) 

where the tilde above the E refers to a average over the 
ostensible distribution. In Eq. (2.15a) we have defined 
a(t — s), this function we label the memory function. On 
a microscopic level it has the form 



a (t ~ s ) = Yl \ gk 



2 e -iO fc (t-a)_ 



(2.16) 



Using the above ostensible distribution we can define 
a linear conditional system state as 



(K}|*(*)) 



(2.17) 



Taking the time derivative and using Eq. (2.7) we get a 
linear differential equation for \ip z (t)} of the from 



dt\j> z (t)) = -r-Hu*(t) + z*(t)L-Ll a(t-s) 



' 8z*(s) 



-ds 



hM*)>, 



(2.18) 



where S/6z*(s) represents a functional derivative. For a 
derivation of this equation see Ref. [10, 16]. The func- 
tional derivative in this equation stops us from calling 



this equation a non-Markovian SSE, as it means that 
c*t|V> z (t)) does not depend upon the state \"ip z (t)) at all 
times for a single function z(t), but rather also upon 
states for other noise functions. That is, we cannot 
stochastically choose z(t) in order to generate a trajec- 
tory independent of other trajectories. Instead, all pos- 
sible trajectories would have to be calculated in parallel, 
which in calculation terms amounts to solving the com- 
plete Schrodingcr equation Eq. (2.7). However, as ex- 
plained in Refs. [11, 12, 16] if we can make the following 
ansatz 

I ^\Mt)) = { %(t,s)\Mt)), (2-19) 
then we can write a linear non-Markovian SSE as 



(2.20) 

where the operator functional (°'F Z (t) is defined as 



(0) iW) 



a(t-s)^f z (t, S )ds. 



(2.21) 



The significance of the superscripts (0) proceeding these 
operators will become apparent in Sec. III. 

To derive the actual (non-linear) non-Markovian SSE 
we need to condition the state on a noise function that 
is equivalent to the actual probability distribution, 

P({a k },t) = (*(*)|{o fc })({o fc }|*(t)). (2.22) 

For most systems ^(i)) is unknown. Nevertheless we 
can use a Girsanov transformation [11, 16] to relate the 
actual noise function to the ostensible noise function. In 
this case, 

z(t) = z A (t)+[ a(t- s)(L) s ds, (2.23) 
Jo 

where z\(t) is equivalent to the noise function used in 
the ostensible case, satisfying the correlations defined in 
Eqs. (2.15a) and (2.15b). With the correct z(t) the actual 
non-Markovian SSE for the normalised state is [11, 16] 



dt\Mt)) 



n 

+((V-(V) t )WF z (t))+z*(t) 



x(L-(L) t ) \Mt)) 



(2.24) 



The notation (L) t is short hand for (ip z (t)\L\i[> z (t)) . From 
Eq. (2.20) and Eq. (2.24) if the operator functional 
(°E Z (t) is known for all time and for each noise function 
z{t) we can solve the coherent non-Markovian SSE. 



2. Quadrature Unravelling- Outlined 

To obtain a non-Markovian SSE with real noise, it is 
natural to consider a quadrature noise operator, 

z{t) = M^'e-^ + St nt (t)e- iw °V*, (2.25) 
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where 6; n t(i) is defined in equation (2.10) and </> is some 
arbitrary phase. The phase <fi defines the measured 
quadrature: an x-quadrature measurement occurs when 
</> is set to zero, and the conjugate measurement of the 
y-quadrature occurs when <fi = k/2. Unless otherwise 
stated we will set 4> to zero. 

The measurement basis for the bath measurement is 
|{<7fe}) and must satisfy 



*(*)!{«*}> = *(*)!{<»}>• 



(2.26) 



The problem with this noise function is that in gen- 
eral it is hard (maybe impossible) to work out a time- 
independent eigenstate in the interaction picture. How- 
ever, we can find this eigenstate if we make the assump- 
tions that for every mode k there exists another mode, 
which we can label —k, such that Q- k = — Q k and 
g-k = g%. These assumptions simply mean that the 
modes coupled to the system come in symmetric pairs 
about the frequency fl. Without loss of generality we 
can take the g k s to be real, absorbing any phases in the 
definitions of the bath operators. With all of these as- 
sumptions we can rewrite equation (2.25) as 

z(t) = ^2 2 9k[X£ cos(Ofet) + Y~ sin(f2 fe i)]. (2.27) 

fc>0 

Here we have introduced the two-mode quadrature oper- 
ators 



X 



k 
I k 



(x k ± X- k )/V2, 
(y k ±y- k )/V2, 

where x k and y k are the quadratures of a k : 

(x k +iy k )/\/2. 



a k 



(2.28a) 
(2.28b) 

(2.29) 



The measurement basis that satisfies Eq. (2.26), in the 
x-quadraturc representation is 



i{%}> = n 

k>0 ' 



r dx' 






J V2^ 


V2 





V2 



„iY, x' 



(2.30) 

With this basis and the above noise operator the noise 
function for the quadrature measurement is 



z(t) 



k>0 



2g k [X+ cos(fi fc i) + Y- sin(Q fe i)], (2.31) 



which by definition is real. 

Furthermore under the above assumptions the memory 
function a(t — s) in Eq. (2.16) reduces to 



(3(t-s) = 2j2\9k\ 2 cos[n k (t- S )}. 



(2.32) 



fc>0 



As in the coherent case we define the ostensible distribu- 
tion as the overlap between the vacuum state and [{<?&}), 
that is 



A({x fc ,y fc }) 



J2k>o( X k +Y k 



1 



(2.33) 



With this distribution the correlation for the real- valued 
noise function is 



E[z(t)z(s)]=P(t-s), 



(2.34) 



where the tilde, like before, means an average over the 
ostensible distribution. For this ostensible distribution 
the differential equation for \ipz{t)} is 



dt\Mt)) = [-~jrH int {t)+z(t)L- L x J (3(t-s) 



5z(s) 



ds 



where L x = L + L'. Making the Ansatz, 



I^W) = (0) 9 Z (M)|^)>, 



(2.35) 



(2.36) 



8z(s) 

the linear non-Markovian SSE is 
d t \Mt)) = l—H^t) + z(t)L - L x ^Q z (t)] \i> z (t)), 



where 



p(t-s)(%(t,s)ds. 



(2.37) 
(2.38) 



To derive the actual non-Markovian SSE we need to 
calculate the correct noise function. The Girsanov trans- 
formation giving the actual real-valued z{t) is [16] 



z{t) = z A (t) + / (L x ) s /3(t - s)ds, (2.39) 
Jo 

where z\(t) satisfies the correlations defined in 
Eq. (2.34). The actual non-Markovian SSE for the 
quadrature unravelling is 



dt|^(*)> 



-H int {t) - (L x - (L x ) t )WQ z (t) 
n 

+ ((L x -(L x ) t )WQ z (t)) t +z(t) 
x(L-(L) t )]\^ z (t)). (2.40) 



Thus, if ^Q z (t) is known for z(t) and all time then we 
can solve the quadrature non-Markovian SSE. 



III. PERTURBATION METHOD 

To solve the non-Markovian SSE, and hence find 
Prod (t) , for the coherent or quadrature unravelling we 
have to work out the operator functionals (°>F z (t) and 
^°'Qz(t) respectively. This has been done exactly only 
for systems for which an analytical solution for p re d{t) 
may be found by other means [11, 12, 14] or for systems 
with a small number of bath modes [16]. In this section 
we to propose our perturbation technique for working out 
these functionals when exact solutions are not possible. 
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A. Perturbation Approach for the Coherent 
Unravelling 

The perturbation that we are going to propose is only 
valid for memory functions of the form 

j 

«(*-*) =X)« W) (*-*)> ( 3 - la ) 

3=1 



where 



v (3') 



(f-a) = IGj- ^e - " 1 * l*--l/a e -*nrf (*-») . (3.1b) 



In principle this is always a valid decomposition for the 
memory function as in the J — >■ 00 and Kj — > limit 
this memory function approaches the microscopic mem- 
ory function displayed in Eq. (2.16). In Rcf. [20] the 
authors suggest that in practice most environments can 
be simulated with J being quite small. 

With this expansion for the memory function the func- 
tional (0) F 2 (i) can be written as 



(3.2) 



< >iW(i) = / a^(t - s) w f z (t, s)ds. (3.3) 



where 



To calculate these operator functionals we set up a set 
of coupled nonlinear differential equations for ^F^ft). 
Taking the time derivative of Eq. (3.3) we get 

d t ^\t) = a«(0) (o >/ z (M)+ ! [d t a^(t-s)j 

Jo 

x Wf M (t,8)ds+ f a U) (t-s) 
Jo 

xd t Wf z (t,s)ds. (3.4) 

The first term is easily evaluated using 

Wf z (t,t)=L, (3.5) 

as derived in Appendix A. The second term is where 
our earlier decomposition of a(t — s) is used. We chose 
a^{t - s) such that d t a^(t - s) oc a^(t - s). This 
results in the second term equaling 



(^+iQ,)<°#«(t). 



(3.6) 



The third term involves the partial derivative 
dt[^fz(t, s)]. To find this we use the fact that 



9, 



5z*(s) 



\Mt)) 



5z*(s) 



dt\Mt)), (3.7) 



which is called the consistency condition in [11]. This 
consistency condition is only valid for t =/= s this is be- 
cause at time t = s the functional derivative is not well 



defined. Using Eq. (2.19) we can write the left-handed 
side (LHS) of the consistency condition as 



5z*(s) 



+ (0) f z (t,s)d t \Mt)). (3.8) 

Substituting Eq. (2.20) in for d t \i> z (t)} gives 
6 



di 



5z*(s) 



\Mt)) 



d t {0) f z (t, S )--Wf z (t,s)H int (t) 
n 

+z*(t)Wf z (t, S )L-^f z (t,s) 
xPm z (t)]\i> z (t)}. (3.9) 



Using Eqs. (2.20) and (2.19) the right-handed side (RHS) 
of the consistency condition gives 



S^ffi^Mt)) = -lH int (t)Wf z (t,s) + z*(t)L 
x(%(t,s)-P^F z (t)(%(t,s) 



^F z (t)]\Mt)). (3.10) 



Sz*(s) 

Equating the LHS with the RHS gives 



d t (0) f z {t lS ) = -^[H int (t),Wf z (t,s)}+z*(t)[L, 



^(t,s)]-[Ln°)F z (t),Wf z (t,s)} 

(3.11) 



-Lt. 



5z*(s) 



{0) F z (t). 



Substituting this equation with Eqs. (3.5) and (3.6) into 
Eq. (3.4) we get 

-[tf<®P z (t),wpv\t)]-V 

k 

(3.12) 

where ^F^' k) (t) is our first order functional. It has the 
form 

^F^' k) (t) = [ a^\t- s) w f { z k \t,s)ds, (3.13) 
Jo 

where we have used the following Ansatz 
S 



5z*(s) 



C°)j7f)(i) = 



(3.14) 



If we knew the form of ^F^' k) (t) then Eq. (3.12) could 
be solved numerically. 
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To find the form of (t) we can take the time The second term as before also simply evaluates to 

derivative of Eq. (3.13). Doing this we get 



dt WP^)(t) = a^(0)W/i fc )(M)+ / [d t a^(t-s)] 

Jo 

x (1) /j fc) (M)ds+ / a u \t-s) 



xd t WfW(t,s)ds 



(3.15) 



The first term is easy to work out. From Eq. (3.12) it 
follows that 



WfW(t,t) = [L ) MFW(t)}. 



(^+<n J )W#«.*)(«). (3.17) 



The third term is worked out via a new consistency con- 
dition, 



dt-^^it) = ^-9 t (°)f>(t). (3.18) 
oz*{s) oz*{s) 



Substituting Eqs. (3.14) and (3.12) into this consistency 
(3.16) condition gives 



9 t W/j fc )(M) = -(y +^ fe )«/i fc )(M) - jr[Hto(.t), ^fi k) (t,s)] + z*(t)[L, S )] 

-[^E*'), {0) ^ k) (t)} - PY (0) ^ w (*). (1) /i fc) (M)i 

_ity^_(i)#(M) (t) . 
^ <5z*(s) z K ' 



(3.19) 



Substituting all these terms into Eq. (3.15) gives 

ftW^)(t) = |G,| 2 [L, ^(i)] - (| + in,) - (f + *O fc ) «^Ci.«(t) 

-|[^nt(*), (1) Fj^)] + z'(t)[L, W^.*)(t)] - p Y {1 WHt), i0) F ( z k Ht)} 

i 

-P Y W F^' k \t, s)} -L^Y {2) F^ k ' l) (t). (3.20) 

l i 

I 



Where the last term is the second order functional, which The differential equation for the n th order functional is 
equals 



^F^\t)= / a U)(t- s )-4-WF(^(t)ds. (3.21) 
Jo Sz s ) 



Here we see that we can develop a general way for 
setting up an n th order differential equations. The n th 
order functional is 



(n)p(j,k,...,l)^ = / a W)( t _ a )(»)/(fc.-.0( t)S )d Sj (3.22) 
Jo 



where we have used the Ansatz 



__(n-l)p(k,...,l) {t) = (n)f(K..J)^ s y ( g 33) 



dt (n)pU,k,...,l)^ = o C3)(0)W/ j (*.-.0( tjt ) 

+ f 1 [ft« (j) (i- S )] ( " , /f "(M)* 



+ / a w (*-*)ft w /i* , '" , ' ; (*»*)d»- 
Jo 

(3.24) 

The first term can always be calculated by the (n — l) th 
differential equation. The second term is always simple to 
calculate as dta^ (t — s) oc (t — s) and the third term 
is always calculable by the (n — 1) order consistency 
condition 



(n) f{k,...,l) i 



Sz*(s) 



Sz*(s) 



(3.25) 



The n th order perturbation method propose is to ter- 
minate this series by setting ( n )Fz' k '"' (t) equal to an 
arbitrary operator. The simplest scheme would be to set 
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this operator to zero, but to keep the theory consistent 
with the Markov limit for all orders, we set 
in the following manner. The zeroth order perturbation 
aries when we use the approximation 



= / a {j) (t- s)dsL. 
Jo 

(3.26) 

Note that the approximation here is the replacement of 
5/5z*(s) by 5/5z*(t). The first order perturbation arises 
when we use the approximation 

ft 

a (j) (t - s)ds[L, {0 F^ k \t)} (3.27) where 



and ^F^\t) is calculated via Eq. (3.12). The n th order 
perturbation arises when we use the approximation 

ft r(n-l)p( fc .— .0 {+\ 
(n)pU, k ,...J ){t) „ / a U) {t _ s)ds d _ W 

= f a ( - j) (t-s)ds[L,^- 1 ^ k '-' l \t)} 
Jo 

(3.28) 

and WpP(t), (»-i)Fp<---< fe )(i) are calculated via Eqs. 
(3.12), (3.20) and (3.24). The physical motivations for 
choosing this type of expansion are; 

a) For most system the memory function will decay and 
thus the most dominant term in the functional derivative 
will be the value as s — > t. 

b) Only ^F^^t) affects the system directly, so the 
further removed the approximation the more accurate 
we expect the approximation to be. 

c) In the Markovian limit, only the zero order term is 
needed. 

To summarize this perturbation method, for environ- 
ments which can be modelled by Eq. (3.1), it is possible 
to obtain a perturbativc solution for the coherent non- 
Markovian SSE. From these SSEs it is possible to gener- 
ate a perturbative solution for /? re< i(i), which by definition 
will always be positive. The number of coupled complex 
differential equations that arc required for this technique 



d (J n +J n +...+J)+d+J = d z J 



J r < 



J -I 



-d+J (3.29) 



where d is the system dimension, J is the number of 
exponentials required to simulate the memory function 
and n is the order of the perturbation. The first term 
represents the number of equations needed to simulate 
the functional derivative. The next term d is for the d 
complex amplitudes of the system. The final term J is 
for the stochastic equations needed to generate the noise 
function z(t). 



B. Perturbation Approach for the Quadrature 
Unravelling 

The perturbation method in the quadrature case is es- 
sentiality the same as the coherent case, but the memory 
function expressed in Eq. (3.1b) is too general. This is 
because the memory function for the quadrature unrav- 
eling must be consistent with the assumptions stated be- 
low Eq. (2.26). The most general memory function that 
satisfies these requirements is 



/3(t- S ) = ]T/^ cos ^- S ), 



(3.30a) 



^ j ' cos \t-s) = 2|Gj| 2 e-^l*- s l/ 2 cos[%(t-s)]. (3.30b) 

This presents a problem as dtP^' COB ^(t — s) is not propor- 
tional to ^' cos \t — s). To get around this we define a 
new function f)w< sin ) (t — s) as 

( t _ B ) = 2|G J fe-^l t - s l/ 2 sin(^(t-s)), (3.31) 

and two functionals 



(0)qO\cos)^ 



ft j - cos \t- s)q z (t,s)ds, (3.32a) 



<°)Q^ sin )(t) = / /3^ sia \t-s)q z {t,s)ds. (3.32b) 
Jo 

The functional ^°'Q z (t) is the found by 

^Q z (t) = J2 {0) Q { z ,cos] ( t )- (3-33) 

3 

Taking the time derivative of Eqs. (3.32a) and (3.32b) we 
get 

dt (°)Q^ cos > (t) = /^'< cos ) (t, t) <°>& (t, t) 

+ ( [d t pV> c ° s \t-s)]Wq z (t,s)ds 



+ / ^ j ' cos \t-s)d t ^q z (t,s)ds, 
Jo 

(3.34a) 



+ / p^-^(t- S )d t ^q z (t, S )ds. 
Jo 

(3.34b) 

As in the coherent case it can be shown that ^°'q z (t, t) — 
L. The two terms involving the derivative of ^' cos '(t—s) 
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and /3w> sm ) (t — s) by definition give 

j d t U ' cos \t-s)^q z (t,s)ds = -^)Q(J^)(t) 

-fi/ 0) Q^' sin) (t) (3.35a) 

f d t p^ (t - s) <■% (t, S )ds = (°)q«-™> (t) 

Jo * 

+^ (0) Q^' cos) (£). (3.35b) 

The last two terms require finding dt^q z (t, s). As in the 
coherent case this is found by the consistency condition 

dtj^-MM) = j^-MMt)), (3.36) 
ozys) oz{s) 

I 



dt (0)ga,cos) W 



where 



(l)Q(,,fc,co S ,cos) W = / 0<j,a»)U a \ O V* W ds> 
Jo ° z ( s ) 

(3.39a) 



(l)gp,sin,cos) W = f'^ift ^^'^ (*) fo. 

Jo M s ) 

(3.39b) 

The higher order functional differential equations are 
found in the same manner as in the coherent case, ex- 
cept the form of [3{t — s) results in 2" as many equations 
for order n. 

The perturbation expansion is similar for this unravel- 
ling, the only difference being that we have 2 n operators 
to approximate. The th order approximation is to set 
the th order functional to 



(o)gjy.«»)( t ) = f ^' cos \t,s)dsL (3.40a) 
Jo 

(o)gO>in)^ = f pU*to)(t, a )daL. (3.40b) 
Jo 



yielding 

d t (0) <z z (M) = -~[H int (t),(%(t,s)]+z(t)[L, 

(% z (t,s))-[LjVQ z (t),Wq z (t,s)] 

-L x -^-^Q z {t). (3-37) 
oz[s) 



Substituting these terms into Eq. (3.34) we get 



I 

The first order approximation is to set the four first order 
functionals to 

(l)gp,COS,CO S )^ = f pU,C°») fa gjfofi, (^Qf.™) (()] , 

Jo 

(3.41a) 

(l)g( 3 ,fc, S i„,co S ) (i) = f p(j,™)fa S } ds [L,(°)Q(>°' C ° S )(t)}, 

Jo 

(3.41b) 

(l)Qb\fc,co S ,sin)^ = f f3 ij ' C ° B) (t,s)ds[L,^Q[ k ' sin \t)}, 

Jo 

(3.41c) 

(DgO-.fe.sin.sin)^ = f pU.^) fa s )ds[L, ^Q^' sin \t)]. 

Jo 

(3.41d) 

and we calculate the th order functionals via Eq. (3.38). 

IV. ENLARGED SYSTEM APPROACH 

To test the accuracy of our perturbation method we 
compare our results for the reduced state with the re- 
duced state found via the enlarged system method of 
Imamoglu [19, 20]. An example of how this method is 



= 2\G 3 \ 2 L- ^p)Q(^-)(t) - ClMQiJ'^fa - ^[H int (t), (°)Q^ cos) (i)] 
+z(t)[L, (°)Q^ cos )(t)] - [L x ^Q z (t), V>Q<t>""\t)] 

-L X J2 (1) Qp< cos ' cos) (t), (3.38a) 

k 

= -^ 0) Q z j ' sin) (t) + %( 0) Qi J ' cos) W ~ -[Hint®, ^Qi 3 ' sin) (t)} + z(t)[L, (°)Qp n )(i)] 

-[Lj 0) Q z {t), m (0) Q^ sin) (0] -L*Y1 (1) QP' sin ' cos) (i), (3.38b) 
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applied to a non-Markovian system can be found in Ref. 
[24]. 

For those who are not familiar with the enlarged sys- 
tem method, we provide a short proof that the reduced 
system dynamics are exactly reproduced by the enlarged 
system method provided that a(t — s), called T(r) in Refs 
[19, 20], is of the form 



«(* - «) = £ \Gj 



2 e ~Hj\t-s\/2-iQ,j(t-s) 



(4.1) 



which is the same as Eq. (3.1). 

The total Hamiltonian for the enlarged system is 



-( 

3 3 

3 

duj \ -kz - %Mct], 

-co V z7r 

(4.2) 

where -ff sys = Hq +H, Cj is the annihilation operator for 
the j th added oscillator and is the Markovian bath 

operator with the correlation 

[0 j (u),i>l(u)] = 6 j , k S( U -u/). (4.3) 

If this is to be the same as Eq. (2.1), then the first two 
lines of Eq. (4.2) must give H sys + -^bath and the final 
line V. Going to the same interaction picture as we did 
in Sec. II A, that is with respect to the Hamiltonians Hq 
and flbath, we get 

VSnt(*) = ihY}G*Le- mt c 3 {rf - Gjtfe^Cjty)]. (4.4) 



Comparing with Eq. (2.9), for the enlarged system 
method to be correct we need 6int (t) = J2j GjCj(t). To 
calculate Cj(t) we use the fact that 

d t Cj(t) = -iujCj(t) - ycj(t) - ^/kJ inj (t) (4.5) 

where v- lri j (t) is the input field which has a time commu- 
tator [Vm,j(t), ^/n ^(s)] = Sj,kS{t — s). For a derivation of 
equation Eq. (4.5) see Ref. [25]. This can be integrated 
to give 



+c,-(0)e 



— Kjt/2 — iujjt 



(4.6) 



It not obvious that J^j GjCj(t) is the same as Eq. (2.10). 
However the time commutator for the bath operators is 



In terms of the enlarged system this means 
£G^[c^)4( S )]e in (<- s ) 

= £ |G-,| 2 e~ Kj ' (t+s)/2 ~ i( ^~ n)(t ~ s) [l 



e + K]( t'+s')/2+iu: 3 (t'~s') 5 ^ _ s /) dt / ds /] 



2 -K,j\t-s\/2—i(u}j-n)(t-s) 



£|G, 



a(t — s), 



(4.8) 



provide a(t — s) has the form depicted in Eq. (4.1). It 
should noted that this result is exact. It is not necessary 
to discard initial transients as in the derivation in Ref 
[20]. 

Since we have shown that the total Hamiltonian for 
the enlarged system is equivalent to the standard non- 
Markovian, then the total states \^Sch(t)) must be the 
same. We can define a reduced state (in the Schrodingcr 
picture) for the enlarged system as Wsch(i) which has the 
Markovian master equation 

dtW Sc h(t) = -~[Hn + H + hJ2uj4 d i + ih J2 

3 3 

x(G*Lc] - Gjfrc^Ws^t)] 

+ £«i©[c,-]Wsch(t). (4.9) 



The reduced state for the system in the fi-interaction 
picture is 

A-cdW = e* Ant Tr en i[W S ch(*)]e-* An * = Tr en i[W Ied {t)}, 

(4.10) 

where the trace is performed over the added oscillators 
and 

(4.11) 

This allows us to define a new master equation for the 
reduced state W ie d(t) as 



-G 3 Dc 3 e-^-^%W, cd {t)] 



(4.12) 



[&int(*),SL(*)K f2(t - s) =a(t-s). 



(4.7) 



which can be solved by standard Markovian techniques, 
for example quantum trajectories [5, 6, 7]. 
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V. NUMERICAL EXAMPLE: THE DRIVEN 
TWO LEVEL ATOM 

In this section we apply our theory to a driven TLA 
with a simple non-Markovian memory function. 

a (t -s) = ^ e iK w -n)(t- s ) e -«|i- s |/2 ) (5 1} 

where tu cnv is the central frequency of the environment, 
k represent the exponential decay of bath memory and 7 
is the Markovian limit decay rate. That is, in the k — > 00 
limit, a{t — s) = "/5(t — s), which is the Markovian limit 
of the memory function [16]. We choose an interaction 
picture such the f2 = u> env so that this memory function 
is simplifies to 

a(t-s) = ^e-^ 2 , (5.2) 

which is consistent with the quadrature unravelings as- 
sumptions. This results in a(t - s) = f3(t — s) . However 
before we apply our theory to the TLA let us revise the 
standard TLA model. 



A. The TLA 

The TLA is one of the most simple quantum systems 
to envisage. It consists of two levels, an excited state 
|e) of energy Tiuj e and a ground state \g) of energy Hlo 3 . 
We define the difference in these energies as huiQ and the 
zero point energy is taken to be the mid point energy 
Ti(uo e + Lu g )/2 = 0. This allows us to define a system 
Hamiltonian as 

H sys = h Y &z ^ 5 ' 3 ^ 

where a z = \e)(e\ — \g)(g\ is one of the spin matrices for 
the TLA. 

Since we are dealing with open quantum systems we 
consider the dynamics of the TLA immersed in the elec- 
tromagnetic field (the bath). In the Schrodingcr picture 
with the dipole and rotating wave approximation (RWA) 
approximation the interaction Hamiltonian is 

V = inJ2(9*k™l-g k d- j a k ), (5.4) 
fc 

where a is the lowering operator for the TLA. This is 
the same form as Eq. (2.3) with L — er, so the above 
non-Markovian SSE theory is applicable to this system. 

If we have a TLA driven by a classical electromagnetic 
field the system Hamiltonian for the TLA under the RWA 
approximation is 

H sys = h^-& z + h^l&e 1 ^ + (5.5) 

where x is the Rabi frequency and Udr is the driving 
frequency of the classical field. However as shown in Eq. 





! 2 4 6 8 10 

time, t [y 1 ] 



FIG. 1: This figure depicts the Bloch vector components of 
the reduced state of a driven TLA calculated by the enlarged 
system method. In this figure all calculations were done using 
the initial system state IV'(O)) = |e) with system parameters 
7 = 1, k=1, x = 5 and A = 3. Time is measured in units 



(2.1) we can also write H sys as Hq + H(t). If Hq = 
Qa z /2, then in the fl interaction picture gives 

H int (t) = n^^a z + h^[ae^ t - nt Ua^e- l ^ t - n % 

(5.6) 

For our purposes we assume = oj^i-- So 

Hi nt (t) = h^&z + f&Vx, (5.7) 

where A = luq — Q is the detuning. 

For the TLA the reduced state can be written in terms 
of the real Bloch vector components x(t), y(t), z{t) as 

Prcd(t) = W + x (t)*x + viQvv + ( 5 - 8 ) 

B. Enlarged System Method 

For the driven TLA with a memory function given by 
Eq. (5.1) the master equation for the enlarged systems is 

d t W Ted {t) = [-—<7 z -—(j x + —{oc-o>c) 1 

W Icd (t)] + K V[c]W Icd (t). (5.9) 

Using 7=1,k = 1,x = 5 and A = 3 the reduced state 
is shown in Fig. 1. For this simple case it was noted 
that the truncation error involved in the enlarged system 
state method was negligible. Because of this we use this 
reduced state for comparison with the ensemble average 
of the non-Markovian SSEs. 



C. Coherent Unravelling-TLA 

Applying the coherent non-Markovian SSE theory to 
the driven TLA, we find that we can rewrite the actual 
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non-Markovian SSE as 
<k\Mt)) = 



* f >t) 

x( )F z (i) + ((at-(at) t )(°^(t); t 
**(t)(o- - <5->t)l IV.C*)), (5-10) 



and the noise function for the TLA becomes 



z(t) = z A (t) + / ait — s){a) s ds. 
Jo 



To calculate the complex amplitudes for the actual 
non-Markovian SSE we apply the system state \ip z (t)) = 
C e (t)\e) + C g (t)\g) to Eq. (5.10) and expand ^F z (t) as 



(5.11) where m = {cr, a\ er 2 , I}. This results 



(5.12) 



d t C g 



d t C e 



l^C g - i^C e + z*C e \C e \ 2 - ^F a , >z C 3 g Cf + ^ z C g \C e \ 2 (l + \C e \ 2 ) 
-^F^ z C 2 C* e {\ + 2\C e \ 2 ) + ^F ItZ C 2 q C* e , 



A .x 

"*yCe — ^ 9 



z*c 2 c: ~ {0 F^ z c e \c g \ 2 {\ + \c e { 2 ) + WF ai . z c 2 c:\c a > 2 



+ m F„^ z C g \C g \ 2 {l + 2\C e \ 2 ) - {0) F Lz C g \C g f 



(5.13a) 



(5.13b) 



In this equation the noise function is given by 

z*(t) = z%(t) + ^e"^ 2 j\^ 2 C g {s)C;{s)d s , (5.14) 

where z A (t) is defined by the correlation 

E[z A (t)z* A (s)} = TL e -^s\l\ (5.15) 

This is generated by having z A (t) obey the following 
stochastic differential equation, 

d t z A (t) = ~z A (t) + ^<:(t), (5.16) 

with z A (0) being a Gaussian random variable (GRV) sat- 
isfying 

E[z A (0)z* A (0)\ = ^ (5.17) 

Here ((t) is standard complex white noise [26] and satis- 
fies E[((t)C(s)} =5{t-s). 

1. th Order Approximation 

For the simple memory function. J — 1, which means 
(°)F 2 (£) = ^F^\t). The th order approximation occurs 
when we assume the form for ^F z (t) in Eq. (3.26). From 
Eq. (5.2) this implies 

^F z {t) = ^{l-e- Kt ' 2 )a, (5.18) 



thus 

(0) F a At) = \{l-er^/ 2 ), (5.19a) 

(0) ^W*) = P ) F <ratZ (t) = <f s >F It ,(t)=0. (5.19b) 
2. 1 st Order Approximation 

The 1 st first order approximation occurs when we as- 
sume a form for ^F^ k \t), by Eqs. (3.27) and (5.2) this 
means 

«F,(t) = 1(1 - e- Kt J 2 )[&, m,(t% (5-20) 

thus 

WF v ,,(t) = 7(l-e- Kt/2 ) (0) i^, z (t), (5.21a) 
(1) F^, z {t) = ^(l-e- Kt/2 ) (0) F a ,At), (5.21b) 
^Jt) = <%/,,(*)= 0. (5.21c) 

The zero order functionals are found by applying the 
TLA operators to Eq. (3.12), giving 

d t m z (t) = l^a-^F z (t) + z*(t)[a,^F z (t)] 

-i[^*, + %**S 0) FM-[* m FM 
^F z (t)]-& ni) F,(t). (5.22) 

Using Eq. (5.12) this gives the following four coupled 
nonlinear equations 



12 



+ (5.23a) 
d t ^At) = -^°'F ff t,(0 + .X^f^W - »A(°»F fft ,W + 2^ !2 (t)[ (0 H,(t) - (0) F ffl , 2 (0] 

- (0) ^t, 2 (t) (0) ^, 2 (t) - - ,*(*)], (5.23b) 
d t {0) F^ z (t) = -^ {0) F*;At) + *f 0) F„iAt) ~ if 0) F*At) ~ {0) F„At)l {0) FiAt) ~ (0) ^W] 

- z *(t)^F a , iZ (t)-^F a , z (t), (5.23c) 

dt^At) = -^ m F Lz (t)-h (1) FaAt)- (5.23d) 



which can be solved in parallel with Eq. (5.13). The zero order functional are given by Eqs. (5.23a) - 

(5.23d). however we now need equations for ( l >F z (t). 
The fist order functionals are found applying TLA op- 
3. 2 nd Order Approximation orators to Eq. (3.20). With a memory function specified 

by Eq. (5.2) we get 

The 2 nd order approximation occurs when we assume a 
form for ^F^' k ' l) (t), by Eqs. (3.28) and (5.2) this means 

~ , dA ] F z {t) = l^[aAF z {t)]~^F z {t) + z*{t) 

< 2 >F,(t) = l{l- e -^)[a^F z {t)l (5-24) 4 

x[a,^F z (t)}-i[-a z + ^a x ,^F z (t)] 

-[A {1 F z (t), {0 F z (t)} - [A 

(\ z (i) = Ttl-e-^V^.W, (5.25a) x^F z (t),^F z (t)] - A^F z (t). (5.26) 

{2) F„,At) = -|(l-e-^ 2 )«F CTt ,,(i), (5.25b) 

(2) ^rt. z (t) = (2) F/, 2 (t) = 0. (5.25c) Using Eq. (5.25) this turns into the four equations 
I 



thus 



d t (1) F a At) = hK (0) F,., z (t) - K^F a , z {t) + lA^'KAt) - *X (1) F. z At) + 2z*(t)V>F amtZ (t) 

+2<% a ,,(*)<% a ,,(t), (5.27a) 
dA ] F ai Jt) = - K WF,i tZ (t) +iA 1) F^ !Z (t)-iA^F a , !Z (t) +2(\, 2 (t)[(°)F Ji2 (i) - (0) ^, z (t)] 

+2^F^ z (t)AFiAt) ~ {1) F„,At)} ~ A ] F^At) ia) F*At) + {0) F a ^ z (t)^F„At)} 
-^F ItZ (t) + ^F a ^ z (t), (5.27b) 

- {0) F^At)} - {0) F„At)l {1) FiAt) - - z*(t) (1) i^t,,(*) 

-| (2) ^(*), (5.27c) 
d* (1) F/, z (i) = -K ( %/,,(t)-| (2) ^,«(*). (5-27d) 



To illustrate how accurate our perturbation method 
is, the difference between the reduced state calculated 
via the enlarged system method and the ensemble aver- 
age from the coherent non-Markovian SSE is plotted in 
Fig. 2. The dotted line corresponds to the th order per- 



turbation, the dashed is the 1 st and the solid is the 2 . 
It is observed that the 1 st and 2 nd order perturbation 
are a lot more accurate then the th order perturbation. 
However, it can be seen that the 2 nd order perturbation 
is not necessarily more accurate than the 1 st order per- 
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2 4 6 8 10 

time, t [y" 1 ] 

FIG. 2: This figure depicts the difference between the re- 
duced state calculated form our perturbative coherent non- 
Markovian SSE and the enlarged system method. The dotted 
line corresponds to the th order perturbation, the dashed is 
the 1 st and the solid is the 2 nd . Other details are as in Fig. 
1. 



turbation. This suggest that our perturbation method is 
an asymptotic expansion rather than a convergent series. 



D. Quadrature Unravelling-TLA 

For the quadrature unravelling the actual non- 
Markovian SSE is 

d t \Mt)) = \ -i^ z -i^a x -{a x -{a x ) t )^Q z {t) 

(o), 



+ ((<7 X - <<7 x ) t ) ™Q z (t) 

+z(t)(a-{a) t )]\^ z (t)), 



and the noise function for the TLA is 



z(t) = z A (t) + f3(t- s){a x ) s ds. 
Jo 



(5.28) 



(5.29) 



Again the coherent case we can calculate the complex 
amplitude equation via applying the state \ip z (t)) = 
C e (t)\e) +C g (t)\g) to Eq. (5.28) and expanding {0) Q z {t) 
as 



(0) Q;W=E™ ( ° )( ^W 



(5.30) 



where m = {a, A, a z , I}. This results in a coupled set 
of differential equations for C e (t) and C g (t) that depend 
on (°^Qm,z(t) and zit). In these equations the real- valued 
noise is given by 



z{t) = z A (i) + ^V Kt/2 / e^ 2 [C g (s)C;( S ) 



where z A (t) is found by 

E[z A (t)z A {s)} 
This is generated by 



7« -«|t-*|/2 

4 



(5.32) 



d t z A (t) = ~z A (t) + |V7CW (5-33) 

with z A (0) being a GRV satisfying E[z A (0)z* A (0)} = 
K7/4. Here £(t) is standard white noise and satisfies 
EmC(s)]=8(t-s) [26]. 



1. th Order Approximation 

The situation is greatly simplified with the mem- 
ory function in Eq. (5.1), as 0(t,s) = f3 (j > cos) (t, s) = 
A j ' cos) {t, s), which in turn implies {0) Q z (t) = 
(°)Q z j ' cos \t) = (°)Qi J ' sin) (t). 

The th order approximation is to set 

C0) 4W = |(l-e- Ki/2 )a, (5.34) 



thus 



(5.35a) 

(0) Q^At) = m Q*„z(t) = i0) QiAt) = 0.(5.35b) 



2. 1 th Order Approximation 
The first order approximation is to set 

(1) i(t) = ^(l-^ /2 )[^ ( °4(t)] (5.36) 

thus 

{1) QaAt) = 7(l-e- Kt/2 ) (0) Q ffz , z (i), (5.37a) 

{1) Q^At) = -^(l-^ Kt/2 ) (0) Q CT t, z (t), (5.37b) 

(1) Q CT t, a (i) - {1) QiAt) = 0. (5.37c) 

The th order functionals are found by applying TLA op- 
erators to Eq. (3.38). With the simple memory function 
this gives 



7K „ K 

■a 

2 



-[A {0) Q z (t), {0) Q z (t)} 

-6r x {1) Q z (t). (5.38) 



-C*Js)C e (s)]ds, 



(5.31) Using Eq. (5.30) this gives, 
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-[ {1) QiAt) + (1) Q^At)l (5.39a) 

dt^Q^At) = -\ (0) ^v(t) + ix {a) Q«;At) -iA(°)Q aV (i) + 2^Q az At)[ {Q) QiAt) - {a) Q^At)} 

- (0) Qat, z (t) (0) Q ff , 2 (t) + (0) Q*t, z W - {1) QiAt) + {1) Q,,At), (5.39b) 
d t {0) Qa^ z (t) = -f (0) Q^W + (0) Q.t, z (*) -f 0) CW - (0) Qa,,W[ (0) g/^W - (0) Q CT ^W] 

+ (0) Q.t, z (t)[(°)Q / ,,(t) + (^, 2 (i)] - zit^Q^A*) 

-\[ (1) Q*At)- (1) Q ff t, a (t)], (5.39c) 

rf t (0) Q^(t) = -~ (0) Q/^W-^ (1) Q^W + (1) Q.t,,W]. (5.39d) 



1l — 

f 0— - 

-1 — 

0.5° 

f 0-— 

-0.5L 

0.5, 
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FIG. 3: This figure depicts the difference between the re- 
duced state calculated form our perturbative quadrature non- 
Markovian SSE and the enlarged system method. The dotted 
line corresponds to the th and the dashed is the 1 st order per- 
turbation. Other details are as in Fig. 1. 



which can be solved in parallel with C e (t) and C g {t). 

To illustrate how accurate our perturbation method is 
for the quadrature unravelling. Fig. 3 shows the dif- 
ference between the reduced state calculated via the en- 
larged system method and the ensemble average from the 
quadrature non- Markovian SSEs for the th (dotted) and 
1 st (dashed) order perturbation. As in the coherent case 
we find the 1 st order perturbation is more accurate then 
the th . 



VI. POST-MARKOVIAN PERTURBATION 

In this section we extend the YDGS post-Markovian 
perturbation [17] to include the quadrature unraveling 
and compare the post-Markovian method with our per- 
turbation method. 

The basis idea behind their perturbation method is to 
expand the operators ^f z (t, s) in powers of (t— s) around 



the point t = s (this is why it is called the post Markovian 
perturbation). That is 

(0) fS,s) = ^Us,s) + [d t ^f z (t,s)\ t=s }(t-s) 

+ \[d 2 t ^f z (t, S )\t=s\(t- S f + ..., (6.1) 

where (°>f z (s,s) = L. To find the first order term we 
simply evaluate Eq. (3.11) at t = s 



d t (0) f z (t, S )\t 



--[H int {s),L]-[Ln°)F z {s),L] 



-D[L, {0 F z (s)]. 



(6.2) 



Thus the functional ^F z (t) for this perturbation is given 
by 

(0) F z (t) = g Q {t)L- gi {t) % -[H- mt {t),L] 

a{t - s)(t - s)[L n °F z {s), L]ds 



a{t- s)(t- s)D[L,^F z (s)]ds, (6.3) 

Jo 

where 



9o(t) 
9i(t) 



a(t — s)ds, 



a(t — s)(t — s)ds. 



(6.4) 
(6.5) 



This equation can not be solved without the initial con- 
dition d t {o) F z (0). However if we make the approximate 
(°F z {s) = Jo a(s - u)Ldu, Eq. (6.3) becomes 

^F z (t) = g (t)L - 9l (t)^[H iat (t),L] - 92(t)[ttL, L], 

(6.6) 



where 
92 (i) 



Jo 



a(t — s)a(s — u){t — s)duds, (6-7) 
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FIG. 4: This figure shows the difference between the reduced 
state calculated from YDGS post-Markovian non-Markovian 
SSE method and the enlarged system method, for both the 
coherent (dotted line) and quadrature (solid line) unraveling. 
Other details are as in Fig. 1. 

which can be solved. The same could be done for the sec- 
ond order terms, but as well as making an approximation 
for (°F z (s) we would need to approximate d s (°F z (s). For 
the purpose of this paper we will only go to first order. 

To extend the idea to the quadrature case we Taylor 
expand the operator ^q z (t, s) in powers of (t — s) around 
the point t = s. To find the first order term we simply 
evaluate Eq. (3.37) at t = s. With the approximation 
(°)Q z (s) = f° 0(s - u)Ldu we get 

{a) Q z (t) = h (t)L - hx(t)~[H int (t),L] - h 2 (t)[L x L,L}. 

(6.8) 

I 



It can be shown that one can establish a set of coupled 
differential equations for these operators provided a(t—s) 
is given by Eq. (3.1). To truncated this perturbation at 
one has to assume a value 0^ n+1 \ It turns out that 
for all operators 0^ n ' other then the only reason 
the operators change from their initial value at t = 
is if the assumed 0^ n+1 ^ is nonzero. This suggest that 
this method is highly dependent on the assumed value 
for 0< n+1 ). 



where 

h (t) = [ /3(t-s)ds, (6.9) 
Jo 

hi(t) = / f3(t- s)(t- s)ds, (6.10) 
Jo 

h 2 {t) = [ [ P(t~ s)/3(s-u)(t- s)duds.(6.U) 

Jo Jo 



For the simple TLA system it is easy to generate these 
approximate expressions for (°'F z (t) and ^°'Q z (t) for all 
time, hence we can obtain solution to the non-Markovian 
SSE. To compare YDGS post-Markovian non-Markovian 
SSE method with our perturbation method, we again plot 
the difference between YDGS method (when 1000 tra- 
jectories where used) and the enlarged systems method. 
The results of this are shown in Fig. 4, where it is ob- 
served that YDGS first order perturbation has a greater 
error than our perturbation method (Figs. 2 and 3). This 
is perhaps not surprising, as the system we modelled has 
k = 1 , which implies it is very non-Markovian. Since one 
of the requirements of YDGS perturbation method is for 
the environment to be close the Markovian regime one 
would expect their method to fail in this regime. 

In Ref. [17] YDGS suggest an alternative perturbation 
method. The functional operator (°'F z (t), which equals 
O z (t) in their notation, is expanded by the functional 
expansion 



(6.12) 



VII. CONCLUSIONS 

In this paper we presented a perturbation method 
for solving the coherent and quadrature non-Markovian 
SSEs. This perturbation method is easily extended to 
any order and is not limited to the post Markovian 
regime. However, the environment is restricted such that 
it has a correlation function satisfying Eq. (3.1). As 
shown in Ref. [20] most non-Markovian environments 
can be simulated via this correlation function with a rela- 
tive small J. This suggest that this perturbation method 
might be useful for simulating non-Markovian evolution 
for p ied (t). 

One appealing feature of this method is that it pro- 



,vi,V2)z(v{)z(v2)dvidv2 + 



o Jo 



(n) (t, v\, ...,v n )z(vi)...z(v n )dv 1 ...dv n + 

I 
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vides a perturbative solution for p ro d(t) which is positive 
by definition. However there is another method, namely 
Imamoglu's enlarged system method [19, 20], which pro- 
vides a better solution for p re d{t)- Imamoglu's enlarged 
system method requires fewer coupled differential equa- 
tions to solve and the only approximation comes in by a 
truncation of the Hilbcrt space of the fictitious modes. As 
one increases the basis size for these modes this method 
will converge to the correct solution. By contrast, con- 
vergence has not been shown for our method. 

This does not mean that our method is useless, as the 
primary interest in our method is not to simulate p Te d(t), 
but to simulate the non-Markovian SSEs. This is in- 
teresting as a continuous in time interpretation of non- 
Markovian SSEs is not clear. In Rcf. [16] we showed 
that these non-Markovian SSE under standard quantum 
measurement theory do not have a continuous measure- 
ment interpretation. However Loubenets in Ref. [22, 23] 
claimed that she has developed a new framework for con- 
tinuous quantum measurements in which non-Markovian 
SSEs represent the evolution of a system state which is 
continuously monitored. 

Future work on this topic is to look into this question. 
Another question that needs answering is whether it is 
possible to derive non-Markovian SSE based on a discrete 
basis such as photon number. We believe this question 
and the previous question will be related. Finally, there 
is the possible application of our method to strongly non- 
Markovian systems such as an atom laser [27] or photon 
emission in a photonic bad-gap material [28, 29]. 



APPENDIX A: DERIVATION OF (0) /*(*,£) = L 

To show that ^°'fz(t,t) — L we start by discretizing 
the functional derivative. We divide the range [0, t) into 



N intervals of width At, so the change in \i[i z (t)) is 



AT-l 



thus 



8z*(s) 



= J2 At 



\Mt)) 



d\MtN)} 



dz*(U)At 



dz*(U), (Al) 



d\MtN)) 
dz*{t t )At : 



(A2) 



if s (ti) is less than t (ijv), which is the only situation 
we are interested in, then taking the limit that s — > t 
(ti = tjv— l) this becomes 

,. S\A(t)) d[|4fev-i))+A^|4fev-i))] , A „, 

lim — -— = ; — ; . (A3) 

s^t 5z*(s) dz*(t N _ 1 )At v ' 

Discretizing Eq. (2.18) we get 



dS z (tN-i)) 



: #mt(<W-l) +Z*{t N - 1 )L 



N-2 



d 



3=0 

x\i> z (t N -i)). 



dz*{tj) 



(A4) 



Substituting this into Eq. (A3) and using the fact that 
the state at time ijv-i only depends on the noise at time 
less then ijv-i, we get the limit 



s\Mt)) 



5z*(t) 

Thus by Eq. (2.19) (0) / z (M) = L. 



L\Mt))- 



(A5) 
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